//
// Created by Administrator on 2021/12/5.
//

#include <cstdint>
#include "MarchingCubes.h"

MarchingCubes::MarchingCubes(RawData *data) : data(data) {
    mesh = new mc_mesh;
    mesh->vertices.reserve(defaultVerticesArraySize);
    mesh->normals.reserve(defaultNormalArraySize);
    mesh->indices.reserve(defaultTriangleArraySize);
}

m_uint MarchingCubes::mc_internal_to_index1D(m_uint i, m_uint j, m_uint k, const mc_vec3i &size) {
    return (k * size.y + j) * size.x + i;
}

m_uint MarchingCubes::mc_internal_to_index1D_slab(m_uint i, m_uint j, m_uint k, const mc_vec3i &size) {
    return size.x * size.y * (k % 2) + j * size.x + i;
}

void
MarchingCubes::mc_internal_compute_edge(mc_vec3i *slab_idx, float va, float vb, int axis, m_uint x, m_uint y, m_uint z,
                                        const mc_vec3i &size) {
    if ((va < 0.0) == (vb < 0.0))
        return;

    mc_vec3f v = {mc_float(x), mc_float(y), mc_float(z)};
    v[axis] += va / (va - vb);
    slab_idx[mc_internal_to_index1D_slab(x, y, z, size)][axis] = m_uint(mesh->vertices.size());
    mesh->vertices.push_back(v);
    mesh->normals.push_back(mc_vec3f({0, 0, 0}));
}

void MarchingCubes::mc_internal_accumulate_normal(m_uint a, m_uint b, m_uint c) {
    mc_vec3f &va = mesh->vertices[a];
    mc_vec3f &vb = mesh->vertices[b];
    mc_vec3f &vc = mesh->vertices[c];
    mc_vec3f ab = va - vb;
    mc_vec3f cb = vc - vb;
    mc_vec3f n = Utils::mc_internal_cross(cb, ab);
    mesh->normals[a] += n;
    mesh->normals[b] += n;
    mesh->normals[c] += n;
}

void MarchingCubes::exec() {
    int nx = data->get_n();
    int ny = data->get_n();
    int nz = data->get_n();
    const mc_vec3i size = {static_cast<m_uint>(nx), static_cast<m_uint>(ny), static_cast<m_uint>(nz)};

    mc_float vs[8] = {0, 0, 0, 0, 0, 0, 0, 0};
    m_uint edge_indices[12];
    auto *slab_idx = new mc_vec3i[nx * ny * 2];

    auto field = data->get_field();

    for (m_uint z = 0; z < nz - 1; z++) {
        for (m_uint y = 0; y < ny - 1; y++) {
            for (m_uint x = 0; x < nx - 1; x++) {
                vs[0] = field[mc_internal_to_index1D(x, y, z, size)];
                vs[1] = field[mc_internal_to_index1D(x + 1, y, z, size)];
                vs[2] = field[mc_internal_to_index1D(x, y + 1, z, size)];
                vs[3] = field[mc_internal_to_index1D(x + 1, y + 1, z, size)];
                vs[4] = field[mc_internal_to_index1D(x, y, z + 1, size)];
                vs[5] = field[mc_internal_to_index1D(x + 1, y, z + 1, size)];
                vs[6] = field[mc_internal_to_index1D(x, y + 1, z + 1, size)];
                vs[7] = field[mc_internal_to_index1D(x + 1, y + 1, z + 1, size)];

                const int config_n =
                        ((vs[0] < 0) << 0) |
                        ((vs[1] < 0) << 1) |
                        ((vs[2] < 0) << 2) |
                        ((vs[3] < 0) << 3) |
                        ((vs[4] < 0) << 4) |
                        ((vs[5] < 0) << 5) |
                        ((vs[6] < 0) << 6) |
                        ((vs[7] < 0) << 7);
                if (config_n == 0 || config_n == 255)
                    continue;

                if (y == 0 && z == 0)
                    mc_internal_compute_edge(slab_idx, vs[0], vs[1], 0, x, y, z, size);
                if (z == 0)
                    mc_internal_compute_edge(slab_idx, vs[2], vs[3], 0, x, y + 1, z, size);
                if (y == 0)
                    mc_internal_compute_edge(slab_idx, vs[4], vs[5], 0, x, y, z + 1, size);
                mc_internal_compute_edge(slab_idx, vs[6], vs[7], 0, x, y + 1, z + 1, size);

                if (x == 0 && z == 0)
                    mc_internal_compute_edge(slab_idx, vs[0], vs[2], 1, x, y, z, size);
                if (z == 0)
                    mc_internal_compute_edge(slab_idx, vs[1], vs[3], 1, x + 1, y, z, size);
                if (x == 0)
                    mc_internal_compute_edge(slab_idx, vs[4], vs[6], 1, x, y, z + 1, size);
                mc_internal_compute_edge(slab_idx, vs[5], vs[7], 1, x + 1, y, z + 1, size);

                if (x == 0 && y == 0)
                    mc_internal_compute_edge(slab_idx, vs[0], vs[4], 2, x, y, z, size);
                if (y == 0)
                    mc_internal_compute_edge(slab_idx, vs[1], vs[5], 2, x + 1, y, z, size);
                if (x == 0)
                    mc_internal_compute_edge(slab_idx, vs[2], vs[6], 2, x, y + 1, z, size);
                mc_internal_compute_edge(slab_idx, vs[3], vs[7], 2, x + 1, y + 1, z, size);

                edge_indices[0] = slab_idx[mc_internal_to_index1D_slab(x, y, z, size)].x;
                edge_indices[1] = slab_idx[mc_internal_to_index1D_slab(x, y + 1, z, size)].x;
                edge_indices[2] = slab_idx[mc_internal_to_index1D_slab(x, y, z + 1, size)].x;
                edge_indices[3] = slab_idx[mc_internal_to_index1D_slab(x, y + 1, z + 1, size)].x;
                edge_indices[4] = slab_idx[mc_internal_to_index1D_slab(x, y, z, size)].y;
                edge_indices[5] = slab_idx[mc_internal_to_index1D_slab(x + 1, y, z, size)].y;
                edge_indices[6] = slab_idx[mc_internal_to_index1D_slab(x, y, z + 1, size)].y;
                edge_indices[7] = slab_idx[mc_internal_to_index1D_slab(x + 1, y, z + 1, size)].y;
                edge_indices[8] = slab_idx[mc_internal_to_index1D_slab(x, y, z, size)].z;
                edge_indices[9] = slab_idx[mc_internal_to_index1D_slab(x + 1, y, z, size)].z;
                edge_indices[10] = slab_idx[mc_internal_to_index1D_slab(x, y + 1, z, size)].z;
                edge_indices[11] = slab_idx[mc_internal_to_index1D_slab(x + 1, y + 1, z, size)].z;

                const uint64_t &config = mc_internal_marching_cube_tris[config_n];
                const size_t n_triangles = config & 0xF;
                const size_t n_indices = n_triangles * 3;
                const size_t &indexBase = mesh->indices.size();
                int offset = 4;
                for (size_t i = 0; i < n_indices; i++) {
                    const int edge = (config >> offset) & 0xF;
                    mesh->indices.push_back(edge_indices[edge]);
                    offset += 4;
                }
                for (size_t i = 0; i < n_triangles; i++) {
                    mc_internal_accumulate_normal(mesh->indices[indexBase + i * 3 + 0],
                                                  mesh->indices[indexBase + i * 3 + 1],
                                                  mesh->indices[indexBase + i * 3 + 2]);
                }
            }
        }
    }
    for (auto &normal: mesh->normals) {
        normal = Utils::mc_internal_normalize(normal);
    }
    delete[] slab_idx;
}
